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We extend the post-processing finite-size (FS) correction method, developed by Kwee, Zhang, and 
Krakauer [Phys. Rev. Lett. 100, 126404 (2008)], to spin polarized systems. The method estimates 
the FS effects in many-body electronic structure calculations of extended systems by a modified 
density functional theory (DFT) calculation, without having to repeat expensive many-body sim- 
ulations. We construct a unified FS DFT exchange-correlation functional for spin unpolarized and 
fully spin polarized systems, under the local density approximation. The results are then interpolated 
to arbitrary spin polarizations. Generalization to other functional forms in DFT are discussed. The 
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sizes demonstrates that it consistently removes most of the FS errors, leading to rapid convergence 
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I. INTRODUCTION 

Many-body simulation methods, such as diffusion 
Monte Carlo (DMC)i~— and auxiliary field Quantum 
Monte Carlo (AFQMC)^ are capable of yielding highly 
accurate results for electronic systems. In extended sys- 
tems, however, because of the use of periodic boundary 
conditions (PBC), finite-size (FS) errors arise in many- 
body (MB) calculations, which are often larger than the 
statistical and other systematic errors. The FS error re- 
flects the spurious interactions of the system, modeled by 
a finite-size simulation cell (supercell), with its own pe- 
riodic images. The long-range nature of the Coulomb in- 
teraction in electronic systems exacerbates the problem, 
and makes the FS error more pronounced. To reduce this 
error, calculations have to be performed using larger and 
larger simulation cells, in order to extrapolate the results 
to the infinite limit of the supercell size. The compu- 
tational cost of MB calculations is typically high. For 
example, QMC methods scale with system size as N 3 - 
N 4 , with large prefactors. Other MB methods typically 
have significantly worse scaling. Thus convergence of MB 
calculations with supercell size is difficult to achieve, and 
accurate extrapolation is often too costly^ FS correction 
schemes that can accelerate this convergence are there- 
fore highly desirable. 

The behavior of the FS error in a MB calculation is dif- 
ferent from that in a standard density-functional theory 
(DFT) calculation^ (or other independent-electron cal- 
culations), as we further discuss in Sec. |TTJ The MB FS 
error consists of a part that is essentially parallel with the 
corresponding FS error in a DFT calculation of the same 
supercell. This part, which will be termed the one-body 
(IB) FS error, is easier to correct. The residual error, 
which will be termed the two-body (2B) FS error, is due 
to the approximate treatment of electron interactions in 
the extended system by the electrons in the simulation 
cell plus their images. The IB and 2B errors are largely 
separable (see Sec.|TT|. The correction of the 2B errors is 



the focus of FS correction methods. 

It has been shown£~— that using a modified periodic 
Coulomb potential instead of the standard Ewald form 
of the electron-electron repulsion can significantly reduce 
FS errors in the calculated ground-state energy and ac- 
celerate convergence. More recently, Chiesa et al^- intro- 
duced a 2B FS correction based on the random phase ap- 
proximation analysis of the momentum distribution and 
2B structure factor, which achieves similar results. Both 
of these approaches are internal correction methods. The 
former requires modification of the form of the electron- 
electron Coulomb interaction, and the latter requires cal- 
culation of the structure factor in the MB calculation. 

The method of Kwee, Zhang, and Krakauer (KZK) 
is an external FS correction approach^ The idea is to 
have independent-electron calculations which mimic the 
FS behavior of MB calculations. KZK developed a FS 
exchange-correlation (XC) functional which describes the 
effect of electron-electron interaction in the spirit of a 
local-density approximation (LDA) in DFT and which 
takes into account the finite size of the simulation cell. 
A functional was parametrized for spin-unpolarized sys- 
tems. The functional has been included in the software 
package Quantum Espresso^ and the method is appli- 
cable to any MB total energy result a 13 ' 15 " — as a simple, 
post-processing correction. 

The FS functional of KZK does not include spin- 
dependence, however. In this paper we generalize 
the method to arbitrarily spin-polarized systems, and 
parametrize a FS XC functional based on the local spin- 
density approximation (LSDA). The new functional in- 
cludes KZK, but will allow FS corrections in systems with 
spin-polarization or magnetic order, for example in many 
transition-metal oxide solids. In such systems the effect 
of electron correlation is generally stronger, and the need 
for MB calculations is greater. The complexity of the sys- 
tems also means our ability to systematically simulate 
larger and larger supercells will be more limited. Ac- 
curate and simple FS corrections will therefore be very 
valuable. In the applications we illustrate how the new 



FS LSDA functional can help accelerate convergence of 
the MB results in various systems. 

The remainder of the paper is organized as follows. 
In Sec. HIl we describe the overall formalism of our FS 
correction and define some notations. In Sec. IIII1 we 
present the FS XC functionals, with the exchange and 
correlation parts each in a separate subsection, and a 
summary of the parameter values of the result of our 
paramctrization in the last subsection. In Scc. lIVl several 
applications are presented to illustrate the FS correction 
method. We further discuss several aspects of the method 
in Sec. |V] before summarizing in Sec. IVII 



II. PRELIMINARIES 

We will consider zero-temperature total energy calcu- 
lations in a cubic supercell of side length L, with M 
atoms in it. Twist boundary condition o 13 ' 21 will usually 
be applied, the many-body generalization of k-points in 
independent-electron calculations. The FS error for a 
given twist k is defined as 

AE{k, L)=E°°- E(k, L), (1) 

where each E denotes per atom energy, and E°° is the 
corresponding result at the infinite supercell limit. 

The FS error can be separated into IB and 2B parts, 
which are also referred to as independent-particle and 
Coulomb FS errors^ The IB FS error arises from in- 
complete k-point integration, i.e., a discrete momentum 
space mesh: 

A£ 1B (k,L) S ]r£(k,L)-i?(k,L), (2) 
M 

where in the first term on the right the sum is over a 
sufficiently large set of k-points (with appropriate weight 
if necessary) to reach convergence. Below we will also 
use the notation AE({k}, L) to denote this term. In a 
standard DFT calculation, because of Bloch's theorem, 
the IB FS error is the only form of FS error: 

#DFT = E BFT (M ^O0,L) = £ D FT({k'} oo,L'), 

(3) 

i.e., AE BFT (k,L) = AE 1B (k,L) in DFT. In Eq. ©, 
different choices of the simulation cell lead to different 
requirements on the set {k} over which to sum, but for 
each choice of L, the infinite limit can be achieved with 
enough k-points. The usual approach is to simply inte- 
grate over a dense k-point grid using the primitive cell. 
The 2B FS error is the residual beyond the IB FS error: 

AE 2B (k,L) ee AE(k,L) - AE 1B (k,L) 

= E°°-E({k},L). (4) 

Thus AE 2B (k,L) is independent of the k twist, i.e., it is 
the FS error aside from k-point sampling errors. 
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FIG. 1: (Color online) Different behaviors of the FS errors in 
DFT and MB calculations, and the separation of the IB and 
2B FS errors. The differences between the calculated AFQMC 
and DFT-LDA total energies in bulk silicon are shown in the 
main graph for three special k-points: V (blue circle), L (red 
triangle), and the Baldereschi22 mean- value k-point (green 
square). The inset shows the calculated ground-state energy. 
The dashed lines are DFT-LDA results, while the solid lines 
are corresponding AFQMC results. The supercell sizes are 
indicated on the horizontal axes. 

This is illustrated in Fig. [TJ which shows the total en- 
ergy difference of bulk silicon between QMC and stan- 
dard DFT LDA calculations. Three special k-points are 
used: T, L, and the Baldereschi 2 ^ mean-value k-point. 
The IB FS errors are seen to be rather parallel between 
the independent-electron and MB calculations. The en- 
ergy difference between QMC and DFT results, the 2B 
FS error in QMC, is essentially independent of k-points, 
consistent with Eq. (gj). The IB FS error in the MB 
calculation can thus be removed by applying a correc- 
tion Ai?Qp T (k, L) from standard DFT. This is frequently 
done in practice. (Note that sometimes the IB FS error 
has a different size dependence compared with that of the 
total AE(k,L), as for the L-point in Fig. [TJ Applying 
this IB FS correction would make the FS error larger^-) 
A residual 2B FS error remains, however, which decays 
slowly with supercell size. 

To consider the 2B FS error in the framework of DFT, 
specifically LSDA, we review how standard LSDA implic- 
itly removes the 2B FS effect, leaving only IB FS error, 
which can be eliminated using sufficiently dense k-point 
grids. The many-electron Hamiltonian in a simulation 
cell is (in Rydberg atomic units): 

H = E v " + E + E yFS ^ - r ^i) . ( 5 ) 

i i i<j 

where i labels an electron, is the electron position, 
the ionic potential on i can be local or non-local, and the 
sum runs over all electrons. We have suppressed the elec- 
tron spin index, but spin is assumed to be present. The 
Coulomb interaction V between electrons depends on 
the supercell size L (and its shape), due to modification 
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by the PBC, for example, via the Ewald method^ The 
corresponding LSDA, as usually formulated, introduces 
a fictitious mean-field system of electrons with the same 
external potential and the same spin configuration, with 
the single-electron Hamiltonianj&ii 

#dft = -V 2 + ^ on (r) + V H {r) + F»(n t (r), n^r), r) , 

(6) 

where the XC potential depends on the electronic den- 
sities for the two spin species, and the Hartree term de- 
pends on the total density. The densities n a {r) are de- 
termined self-consistently, with Fermi statistics imposed 
on the occupied orbitals (eigenfunctions of ^dft)- The 
LSDA XC functional in Eq. has an oo in the super- 
script: V££ is supposed to be the XC functional of the 
infinite system. 

The key is that using V££ removes any 2B FS effect. 
Typically is obtained from QMC results on the ho- 
mogeneous electron gas (HEG), extrapolated to infinite 
size^ ' 23 ' 24 Because LSDA yields an independent-electron 
problem, in which the Hamiltonian is only dependent on 
the electronic densities, Bloch's theorem can be applied. 
The only requirement for size convergence is the con- 
vergence of the densities, which are fully defined by the 
periodicity of the system once the primitive cell is spec- 
ified. The extrapolation of FS HEG data to the infinite 
size was important in that it eliminated the size depen- 
dence in the XC functional, making the theory consistent 
at the thermodynamic limit. 

To have an LSDA theory for the FS supercell, we must 
recover the FS dependence of the interaction term, and 
make it parallel to the MB Hamiltonian: 

ffDFT = -V 2 + V ion (r) + V h (t) + K F c S K(r), H(r), r) , 

(7) 

where the FS V^, s represents the XC effects in the finite 
simulation cell in which the MB calculation is carried out. 
Within DFT, the XC energy functional e^ c is constructed 
based on the calculations of the HEG, and then applied 
to realistic system by the LSDA approximation* 1 * 2 ^— 
Below, we examine the electron gas in finite cubic super- 
cells (of linear size L) to obtain an ansatz for the FS XC 
energy e xc (L; n-t-, nj.), from which we can derive the FS 
XC potential V^ s . A functional will be parametrized as 
a function of L. With such a functional, LSDA calcula- 
tions can be carried out in the same supercell as the MB 
calculation. The difference between the FS LSDA and 
the infinite LSDA results (e.g., from usual dense k-point 
calculations), A£'p| T (k, L), will provide an estimate of 
the MB FS error. The FS corrections used in the remain- 
der of the paper can be summarized as: 

ADFT 1B = £ DFT (oo) - E DFT (L) , 
ADFT 2B = E DFT (L) — E™ T (L) , (8) 
ADFT FS = ADFT 1B + ADFT 2B . 



III. FS XC FUNCTIONAL 

In this section we discuss the FS XC functional in the 
HEG. We will use N to denote the total number of elec- 
trons in the supercell of volume f2 = L 3 , and N-f and iVj. 
for spin-f and spin- J,, respectively. For convenience, we 
will refer to the majority spin as t _ spin below, i.e., as- 
suming JVf > N^. The density is as usual specified by r s , 
the radius of a sphere containing one electron on average, 
which is related to the charge density by 47rr 3 /3 = l/n. 
The fractional spin polarization parameter £ is defined 
as the ratio of spin density and charge density: 



where n-j- (njj is the electronic density of spin compo- 
nent up (down). The total charge density is n=n^+n^. 
The spin unpolarizcd state corresponds to £ = 0, while 
the fully spin polarized state to £ = 1. The desired FS 
XC functional e xc (L; n-)-, njj is thus equivalcntly writ- 
ten as e xc (r s , L, () in the HEG. The total XC energy 
is e xc {r s ,L,C) = £ x {r 8 ,L,Q + £ c (r s , L,(). We discuss 
the exchange and correlation energies separately in See's 
IIII Al and IIIIBI and summarize our parametrized func- 
tional with numerical parameter values in Sec. IIII CI 



A. Exchange energy functional 

The exact total energy of the HEG can be separated 
into the HF energy, which is the sum of kinetic energy 
and the exchange energy, and the remainder, which is 
termed correlation energy. The kinetic energy is simply 
given by a sum of the energies of independent electrons 
filled from the lowest single-particle energy state up to 
the Fermi surface. The exchange energy from the filled 
Fermi sphere can be calculated analytically, and is of the 
form cx l/r s . The asymptotic expression of the corre- 
lation energy in the high density region was derived by 
Gell-Mann and Brueckner^i which is used in combina- 
tion with DMC results at various densities^ to obtain 
a parametrized fit 2 ^ of the LSDA XC energy As men- 
tioned, all of these are at the infinite size limit . 

In finite-sized supercells, the HF exchange energy of 
the HEG scales as 1/L. The FS error in the exchange 
energy per electron, Ae x (r s ,L) = e£°(r s ) - e x (r s ,L), 
scales^ as l/r s , with a coefficient Av(N) which is only 
dependent on N (i.e., the ratio of L/r s ), not on r s . The 
FS exchange energy approaches the infinite size limit 
from below. Av(N) decays smoothly like l/N 2 / 3 when 
averaged over many k-points^ 3 - Thus the first order cor- 
rection of exchange energy will be of order r s /L 2 . We will 
include a next-order term in the parametrization, r 2 /L 3 , 
which follows the same scaling relation. 

In a system with two spin components, each compo- 
nent contributes independently to the total exchange en- 
ergy There is no exchange interaction between opposite 
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spins. Thus, for a system of polarization £, the exchange 
energy is given by 

e*(r„C) = ^£x(r st ,0) + i-^e x (r si , 0) (10) 



where £:r(r s , 0) is the exchange energy of the unpolarizcd 
system with total density given by r s , and the density for 
the spin up (down) component is given by r s f = r s (l + 
C)" 1/3 (r si = r s (l - C)~ 1/3 )- With this relation, we 
can obtain the FS exchange energy for an arbitrary spin 
polarization. 

However, the FS exchange energy for spin unpolarizcd 
system is discontinuous as given by KZK.— At low den- 
sities, the exchange energy is forced to go to zero rapidly, 
because in the finite supercell there will only be less than 
one electron in each spin component when r s is larger 
than some threshold j x . In this regime, the functional 
form from scaling, i.e., r" _1 /i™, would diverge. A dif- 
ferent functional form is used instead, and the two are 
joined at r s = j x by requiring the exchange potential to 
be continuous, which resulted in a discontinuity in the 
exchange energy. With this form of the exchange energy, 
the use of Eq. (|T0|) would lead to multiple discontinuities 
in a partially polarized system, and an exchange that 
tends to be too high in the intermediate densities which 
connect the two density regimes. 

Instead of the formula in Eq. (jTU)) , we will make an in- 
terpolation between spin unpolarizcd and fully polarized 



systems. The interpolation formul; 
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e x (r a , L, = e x (r s ,L, 0)+f(Q[e x {r s ,L, l)-e x (r a , L, 0)], 

(11) 

with the function /(£) defined as 

/(o= (i+c)4/ ;:, (i _- 2 c)4/3 - 2 - c» 

The formula in Eq. (jlip gives the exchange energy for any 
polarization by interpolation bwteeen the exchange ener- 
gies for unpolarized and fully polarized systems. When 
L is infinite, this formula reduces to the exact relation in 
Eq. (|T0|) . The use of this formula also makes the treat- 
ment of the exchange consistent with that of the corre- 
lation (see Sec. MI B[) . for which there is no analogy of 
Eq. CEUD. 

We will parametrize the exchange energies for the two 
end point systems in the following form: 



a o<J>) _|_ a i(p) „ j_ Q 2(p) „2 



e x (r s ,L,p) 



a.3{p)L , q 4 (p)L . a 5 (p)L' 



other. 



(13) 

Here p has two discreet values: p = for spin unpolarized 
while p = 1 for fully polarized. The formula is the same 
as that in Ref. [TH for the high density region, but we have 
introduced two additional terms in the low density part. 
The values of ao(p) are determined by the infinite size 



HEG exchange energies. The other coefficients arc ob- 
tained from a fit and by continuity conditions at 7 X . As 
described in Sec. MI C| wc choose the same j x for p = 
and p = 1, hence for all polarizations £. At r s = y x , we 
require the exchange potential and its first derivative to 
be continuous for all polarizations. The conditions are 
straightforward to impose at p = and p = 1. For an 
arbitrary spin polarization, the spin interpolation func- 
tion /(£) in Eq. (jTTJ) depends on the spin density. The 
exchange potential is given by 



d(ns x (r s ,L,()) 



dn ni) 

v x (r s ,L,0) + f(()[Mrs,L,l) - P x (r s ,L,0)} 

„ 3/(0 



dn 



t(4) 



-[e x (r a ,L, 1) - e x (r s ,L,0)]. 



(14) 



We see that, if the energy difference between the polar- 
ized and unpolarized systems, s x (r a , L, 1)-E x (r s , L, 0), is 
made continuous at r s = y x , the exchange potential will 
be continuous for any It is easy to show that conti- 
nuity of the first derivative of the exchange potential is 
also satisfied with these conditions. We calculate the ex- 
change energies for different supercell sizes by averaging 
over k-points. The quality of the fits is consistent with 
that in Ref. [TH Numerical values of the parameters are 
given in Sec. MI CI 



B. Correlation energy functional 

Unlike the exchange energy, the correlation energy con- 
sists of both intra- and inter-spin contributions. There 
is no exact mapping of the correlation energy at an arbi- 
trary polarization from spin unpolarizcd results. Within 
standard DFT, the correlation energy for partial spin po- 
larization is obtained from an interpolation between spin 
unpolarized and fully spin polarized systems, as we have 
done with the exchange energy. We will follow the same 
procedure. Below, we first construct the energy functions 
for the two end-point systems of p = and p = 1. 

In Ref. [T3I the correlation energy of unpolarized system 
of the HEG is expressed separately for high, intermedi- 
ate, and low density regions. In the high density region, 
the formula is parameterized by fitting to Ceperley and 
Alder's extrapolation fiti^"— 



e c {r s ,L)=e™(r s )-Ae x (r s ,L) 

+ [h(r s ) - l]Ae K (r s ,L) + b 2 (r s )/N, 



(15) 



where the parameters b\ (r s ) and 62 (r s ) were written as 

^ /2 1/2 

low order polynomials in r s and r s , respectively, 
and fitted to tabulated data for unpolarized systems. 30,31 
(Note that here Ae x (r s ,L) and Ae,R-(r s ,L) are given by 
£ 2; (r s ,L)-£^ (r s ) and e K {r a , L)-e<g(r s ), respectively.) To 
avoid the divergence at sufficiently large r s or small N, 
the correlation energy is set to zero. In the intermedi- 
ate density region, a polynomial function is used whose 
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FIG. 2: (Color online) Quality of the new FS correlation func- 
tional, a) The new parametrization of the correlation energy 
functional for the spin unpolarized HEG is compared with 
AFQMC data and with the previous functional of KZK.— 
The solid lines are the new form given by Eq. I)16[l. while 
the dashed lines are from KZK. b) The parametrized correla- 
tion energy for fully spin polarized HEG systems is compared 
with AFQMC data. Solid lines represent the fits while sym- 
bols are from AFQMC at finite sizes. Statistical error bars 
in the AFQMC results are smaller than the symbol size. The 
infinite size LSD A result is also shown for p = 1. 



parameters were completely determined from continuity 
conditions. For fully polarized systems, the available 
finite-size data from QMC were not as detailed. We 
instead modify the procedure above, by generating FS 
HEG data in the intermediate density regime and fitting 
them to obtain a unified function extending to higher 
densities. 

We parameterize the spin unpolarized and fully spin 
polarized correlation energies by a new unified function 
with only two density regions: 



'(r a ,p) 



~L,p) = 



a i(p) , 

L 2 1 



g(r s ,L,p) 
L 3 



0. 



if r s < 7 C 



other. 



(16) 

The functional form is chosen to approach the correct 
asymptotic value at high density, where we have used the 
Perdew-Zunger (PZ) parametrization^! for the function 
e^(r s ,p). The leading FS term ai(p) exactly cancels its 
counterpart in the exchange energy to ensure that the 
total XC energy scales as 0(1/L 3 ) as expected. For the 
1/i 3 term, we choose the functional form 



g(r s ,L,p) = g 1 (L,p)r s hi(r s ) +g 2 (L,p)r s + g 3 (p)rl /2 
+ 9i{v) r 3 J 2 In(r s ) + g 5 (p) r z J 2 + g 6 (p) r 2 s . 

(17) 

The last four parameters arc determined by the fits, while 
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FIG. 3: (Color online) Interpolations of the FS correlation 
energy for partial spin polarizations in the HEG. Three su- 
percell sizes are shown, each with two lines showing the two 
interpolation schemes. Corresponding FS AFQMC results are 
shown together with statistical error bars. 



the other two are determined by continuity conditions. 
Wc require the correlation energy and its first derivative 
to be continuous at J c (p)- Our choice of j c (p) and the 
final parameter values are given in Sec. MI CI 

We perform phaseless AFQMC^ calculations of the cor- 
relation energies of the HEG, for p = 1 and p = 0, in 
finite supcrcclls at typical densities. With each choice 
of N and L, we average the total energy over multiple 
k-points. The results are shown in Fig. [2l it is seen 
that the new parametrized functions fit all the AFQMC 
data well. Figure [2^l also shows a comparison between 
the new functional and that of KZK for unpolarized sys- 
tems. The two fits are almost indistinguishable, except in 
the intermediate region, where the new function shows a 
slightly better fit to AFQMC results, especially for small 
lattice sizes. As mentioned, in the approach of KZK, 
the functional form at high density is fully determined 
by existing DMC data. AFQMC data were only used to 
guide the choice of the boundaries which separated the 
high density regime from the intermediate and low den- 
sity regimes. In contrast, the current approach only had 
AFQMC data at intermediate densities to determine the 
functional form for r s < j c . The good agreement between 
the two paramctrizations across the density range is thus 
very reassuring. Contrary to the exchange energy, the 
correlation energy approaches the thermodynamic limit 
from above, as illustrated in Fig. [2b . 

We have now obtained the FS correlation energy for- 
mulas for spin unpolarized (p = 0) and fully spin polar- 
ized (p = 1) systems, with which we can interpolate to 
arbitrary spin polarization £. As mentioned, the interpo- 
lation of the correlation energy is not exact. In standard 
DFT, two forms have been widely adopted) ^ 26 ' 32 that 
of PZ as we have used for the exchange in Eq. (jTTj) . and 
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TABLE I: Coefficients in the parametrized FS LSDA XC func- 
tional for spin unpolarized and fully polarized systems. All 
values are given in Rydberg atomic units. The choices of the 
continuity densities 7^ and 7 C are given in Sec. IIII CI 
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p=0 


V- 




ai(0) 




di{l) 


5.(1) 





-0.9163 




-1.1545 




1 


-2.2037 




-1.7491 




2 


0.4710 




0.2967 




3 


0.2339 


0.2109 


0.1812 


0.7528 


4 


-0.4880 


8.4987 


-0.4515 


3.3314 


5 


0.1847 


-13.6840 


0.1786 


-5.1050 


6 




-4.6977 




-2.3048 



that of Perdew and Wang (PW): 

s c (r s , L, C) = e c (r s ,L, 0) + a c (r s )-^[l - C 4 ] 

+ [e c (r s ,L,l)-e c (r Sl L : 0)}f(C)C\ (18) 

where a c (r s ) is the spin stiffness function which in- 
cludes additional fitting parameters. The interpola- 
tions for several different supcrcell sizes are shown in 
Fig. [3l In the PW results, we have fixed the a c (r s ) 
at the infinite-size values without further parametriza- 
tion. It is seen that the two schemes give similar re- 
sults. Compared to the actual FS data from AFQMC 
calculations in polarized systems, both schemes provide 
reasonable estimates. Discrepancies are visible, espe- 
cially at smaller supcrcell sizes. In the applications 
below, we have used the PZ interpolation scheme for 
the correlation energy, due to its simplicity, i.e., us- 
ing the same interpolation as in the exchange energy in 
Eq. (fTTj) . The PZ form is equivalent to the special choice 
a c (r s ) = f" (0)[e c (r s ,L,l) - e c (r s ,L,0)} inEq.[]l This 
is further discussed in Sec. |V] 



C. Numerical Parameters 

The functional forms of our parametrized FS LSDA 
XC functional are given in Eqs. (H3J), (HU), and ([T7]). By 
performing AFQMC calculations and fitting to the re- 
sults with the procedures described above, we obtain the 
values of the coefficients, which are given in Table |U The 
only remaining parameters are the densities at which the 
different pieces in the functional are joined. For the ex- 
change energy, we used 7^ — I'siN = 1) for both spin un- 
polarized and fully polarized states. For the correlation 
energy, we chose j c = r s (N = 1) for the fully polarized 
system and j c = r s (N = 0.5) for the unpolarized system. 

These choices are guided by the basic idea that beyond 
7, we have a supercell of the HEG in which there is only 
a fraction of an electron. The exchange (correlation) en- 
ergy needs to vanish rapidly in this regime. Clearly the 
location of 7 cannot be uniquely defined. The results 
from the XC functional should not be sensitive to the 



precise value of 7. We have verified that this is indeed 
the case in typical supercell sizes. In principle, j x should 
depend on the spin polarization For the fully polar- 
ized system, N = 1 signals the transition point, while for 
the unpolarized system, N = 2. We have taken them to 
be the same value, which simplifies the interpolation and 
makes the functional easier to implement, because all val- 
ues of C share the same j x . The form of the correlation 
energy is continuous and it is easier to have different 7 
values for different p. We have thus chosen 7 C for p = 
to be the same as in KZK for consistency. 

In the correlation energy, there are six parameters in 
the formula of Eq. (fTTj) . We only need to fit the last four, 
since the other two, 171 and g2, are fully determined by 
the continuity conditions. This is why they have a de- 
pendence on lattice size L. In the applications below, we 
have used the PZ interpolation scheme in the correlation 
energy, i.e., using the same interpolation as in the ex- 
change energy. The 1/L 2 term is then exactly canceled 
at any polarization, leaving only l/f2 and higher order 
terms in the FS errors. 



IV. APPLICATIONS 

We consider four examples here to illustrate the appli- 
cation of the new FS LSDA functional. The four prob- 
lems are chosen to cover different aspects of the appli- 
cation. All of them have spin polarizations or magnetic 
ground states which require the new spin-dependent FS 
functional. Calculations in single atoms and molecules 
provide systematic MB results for arbitrary supercell 
sizes which allow a detailed study of the FS convergence. 
The ionization energy of Mn is calculated to study the ef- 
fect of the FS LSDA in supercells when there is a physical 
charge imbalance. In the last example, transition metal 
oxide systems MnO and FeO are considered to study the 
application of the FS correction in extended systems with 
magnetic order. In the first three examples, we use the 
phaseless AFQMC method^ 3 - 3 - in a plane- wave basis, with 
a norm-conserving pseudopotcntial to carry out the MB 
calculations, before applying the FS corrections. More 
systematic AFQMC calculations of atomic and molec- 
ular systems using plane-wave basis exist^^ Our goal 
here is simply to use it as an accurate MB approach for 
model extended systems to study FS corrections. In the 
last example, we use existing DMC result o 34 ' 35 and use 
the new FS functional as a post-processing correction. 

Our first application is to calculate the dissocia- 
tion energy (D e ) of the P2 molecule, using a plane- 
wave basis and PBC with single k-point sampling, 
k = T. In both AFQMC and DFT calculations, a 
norm-conserving Kleinman-Bylander— separable nonlo- 
cal LDA pscudopotential 3 -^ is used. The P atomic 
ground state is magnetic, with an electronic configura- 
tion 3s [t|] 3p[ttt]i a "close-shell" system for each spin. 
The P2 molecular ground state is non-magnetic, with 5 f- 
spin and 5 J,-spin electrons. The total energy calculations 
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FIG. 4: (Color online) Convergence of the P2 molecular dis- 
sociation energy vs. supercell size. Calculated AFQMC total 
energies of the P atom and P2 molecule are shown in pan- 
els (a) and (b) , respectively. The corresponding D e of P2 is 
shown in (c). The infinite-size limits of DFT and AFQMC 
are shown by the dashed and dot-dashed lines, respectively. 
AFQMC statistical error bars are shown but many are smaller 
than the symbol size. The orange arrow in (c) points to the 
experimental value. 



are performed in cubic supercells of size L = 7— 18 Bohr. 
Figure Q] shows the MB results from AFQMC calcula- 
tions, as well as DFT results from both the standard, 
infinite-size LSDA XC functional (denoted by DFT°°) 
and the new FS functional (denoted by DFT FS ). 

The conventional DFT°° energies converge to the in- 
finite limit very rapidly, as expected. The DFT FS cal- 
culations exhibit a significantly slower convergence, re- 
flecting the 2B FS effect which it is designed to capture. 
For the P atom, the overall FS effect is not very large, 
but a persistent FS error can be seen in the raw AFQMC 
results. Up to L = 18 Bohr, the energy is still about 
0.15 eV away from the infinite-size value. The IB cor- 
rection with standard DFT [ADFT 1B ; see Eq. ©] is in 
the wrong direction and leads to a larger FS error and 
slower rate of convergence. By contrast, the DFT FS to- 
tal energies, which include both IB and 2B corrections 
[Eq. ((HJ], show size dependence similar to the AFQMC 
results. After ADFT FS correction, the AFQMC total en- 
ergy converges to the asymptotic value rapidly, as seen 
from the essentially flat curve beyond L = 9 Bohr. 

A similar trend is seen in the P2 molecule as shown in 
Fig.HjD. In this case the IB correction does reduce the FS 
error, but only up to intermediate size L, beyond which 
the standard DFT has reached convergence. The MB 
results remain unchanged for large L, where the 2B FS 
effect dominates. The DFT FS calculation again tracks 
the FS dependence of the MB AFQMC, and the total 



FIG. 5: (Color online) Convergence of the Si2 molecular dis- 
sociation energy vs. supercell size. Same conventions as in 

Fig.rj 



energy of the latter converges very rapidly to the infi- 
nite size limit after applying the full FS correction. The 
P2 system is unpolarized, so the KZK parametrized FS 
functional can also be applied. Compared to the results 
there^ the new functional seems to slightly overestimate 
the FS error. This small difference is from the choice 
of exchange energy discontinuity point j x , as mentioned 
in Sec. IIII CI With the FS corrections, the calculated 
D e shows good convergence by about L ~ 12 Bohr, as 
seen in Fig. Hfc. 

Our second application is for the Si atom and molecule. 
These are "open-shell" systems, and both systems are 
spin polarized in the ground state (3 j 1 I and 5 j 3 |, 
respectively). The calculations are done in a similar way 
to the first application, with planewave basis and a norm- 
conserving non-local pseudopotential in a periodic super- 
cell. The dependence of total energies and the dissoci- 
ation energy with lattice size L is presented in Fig. [5] 
Similar to the P atom, the usual DFT IB correction for 
the Si atom has an opposite sign from the total FS er- 
ror. With DFT FS corrections, the convergence of the 
total energy of the atom is greatly improved across the 
entire range of lattice size studied. Similarly, for the Si2 
molecule, the total energy becomes essentially fiat when 
the full correction is included, leading to rapid conver- 
gence to the infinite-size limit. 

Although the atoms and molecules are relatively sim- 
ple systems for MB calculations, they provide excellent 
tests for the FS correction method. They are model sys- 
tems of a molecular solid in which the supercell size can 
be arbitrarily and systematically varied. Because of the 
nature of the interactions and the highly inhomogeneous 
density distributions in the supercell, LSDA calculations 
are not particularly accurate, as can be seen by the final 
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DFT, Makov and Payne (MP)i^ argued that the de- 
pendence has the form: 
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FIG. 6: (Color online) Convergence of the ionization energies 
of the Mn atom vs. supercell size (plotted as the inverse of 
the supercell volume 1/n). Panel (a) shows the calculated 
AFQMC total energies for Mn, Mn + , and Mn ++ ; the corre- 
sponding first (IP) and second (IIP) ionization energies are 
shown in panels (b) and (c). Curves labeled QMC are from 
raw AFQMC results [for the charged ions, the leading order 
AMP (1) correction is included (see text)]. The ADFT 1B and 
ADFT FS corrections are as in Figs. [4] and \5\ [for the charged 
ions, these corrections also add the AMP' 2 ' correction (see 
text).] Statistical error bars are shown and are about the size 
of the symbols. The heavy black dash lines indicate the infi- 
nite size limits of the calculated results in panel (a) and the 
experimental ionization energies in panels (b) and (c). 



DFT D e 's (panel (c) in Figs. Hand 0}. Yet the FS correc- 
tion method based on LSDA works extremely well. This 
underscores an important point, namely that for the FS 
correction method to be effective, the system does not 
necessarily have to be weakly or moderately correlated. 
The FS correction method requires LSDA to provide a 
good approximation in capturing the difference between 
the system with interaction V FS and that with the full 
interaction V (no PBC image interactions)* 1 ^ This is not 
the same as requiring LSDA to work well in either system. 
As the supercell size increases, the correction decreases 
and approaches zero. One would expect the correction 
to be ineffective if the supercell is less than the size of 
the XC hole, which would cause a distortion that is not 
captured in any LSDA approximation. 

The third application is to the Mn transition metal 
atom and its first and second ionization potentials. The 
k-point used here is (0.25,0.25,0.25). For charged sys- 
tems such as Mn + or Mn ++ , there is an additional com- 
plication in a supercell calculation, since charged systems 
are incompatible with PBC. A uniform neutralizing back- 
ground is introduced to maintain charge neutrality— but 
this results in a slow 0(1/L) size convergence to leading 
order i 39 ' 40 In independent-electron calculations such as 



E(L) = E(oo) 



2 

q a 
~2L 



2nqQ 



0(ir 



(19) 



where q is the neutralizing charge, a is the Madelung 
constant, and Q is the quadrupole moment of the su- 
percell. Thus, adding the counter-terms AMP^ = %^ 
and AMP( 2 ) = -^gs 2 to E(L) accelerates the size- 
convergence of DFT calculations for charged systems. 
The AMP' 1 ' correction applies equally well to a many- 
body calculation with PBC, since it depends only on the 
neutralizing charge, and all AFQMC results reported be- 
low for charged systems include this important correc- 
tion. The AMP< 2 ) correction depends on Q, which will 
generally differ, in a many-body calculation, from the 
DFT value. In this work we have not calculated Q in 
AFQMC, but have used the DFT value. 

Figure [5] shows the size convergence of total energies 
of Mn, Mn + , and Mn ++ and the corresponding first and 
second ionization energies. For Mn and Mn + total ener- 
gies, ADFT FS removes nearly all of the FS error, while it 
somewhat overcorrects for Mn ++ . The ADFT 1B correc- 
tions are smaller in all cases. For the ionization energies, 
both ADFT 1B and ADFT FS are seen to overcorrect the 
FS error. For IP, the residual error after ADFT FS cor- 
rection is slightly smaller than the raw AFQMC result, 
while ADFT 1B increases the FS error. For IIP, both FS 
corrections increase the FS error. Extrapolation of these 
results to the infinite size limit, leads to IP and IIP val- 
ues of 7.27(08) and 23.15(08) eV, respectively. These are 
in reasonably good agreement with the experimental val- 
ues of 7.43 eV and 23.07 eV, respectively^ Because of 
strong FS error cancellations, the raw QMC IP and es- 
pecially IIP energies show little FS effect. This magnifies 
the residual errors in the FS correction for charged sys- 
tems. The most likely source of the residual error would 
appear to be the discrepancy in applying the AMP^ 2 ' 
directly from DFT. The MB and DFT estimates of the 
quadrupole moment Q are likely to differ. We will leave 
further investigation of this effect to a future study. The 
total energies, after FS correction, show much smaller 
FS error, especially in the neutral system of Mn. Given 
the strong magnetic nature of this transition metal sys- 
tem, this is an encouraging result for our FS correction 
scheme. 

The last application is for crystalline MnO and FeO. 
These materials are challenging to simulate due to 
the complex interplay of strong correlation effects, p- 
d hybridization, and magnetism. We apply the DFT 
FS correction directly to previous results from DMC 
calculations! 34 ' 35 Although the pseudopotentials used in 
our modified DFT calculations arc different from those 
used in DMC, this has only a weak effect on the FS 
corrections. As seen in Fig. [3 the ADFT FS correction 
greatly reduces the DMC FS error. In both materi- 
als, it somewhat over-corrects and approaches the infi- 
nite limit from above. The structure factor correction of 
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FIG. 7: (Color online) The relative FS errors (per number of 
formula units Z in the supercell) for a) MnO in the Bl phase, 

° 3 

and b) FeO in the iB8 phase at volumes of 21.7 A /MnO and 

° 3 

20.4 A /FeO, respectively. Results from the DMC calcula- 
tions are indicated by green circles; corrected values using the 
structure factor method (taken from Ref. [34ll35h are shown 
by the blue squares. DMC results corrected by ADFT FS are 
indicated by red triangles. 

Chies a 12 ' 34 ' 35 , also shown in tho figure, somewhat under- 
corrects. For MnO the residual error is somewhat larger 
than for ADFT FS , while it is smaller for FeO. A possi- 
ble difficulty for ADFT FS is the irregular supercells used 
for these materials. The ADFT FS correction method is 
based on cubic supercells, parametrized by a single length 
scale L. This is a good approximation for other high 
symmetry systems such as fee supercells.— If the super- 
cell shape is far from cubic and has strong anisotropy, 
however, the FS total energy will increase, resulting in a 
smaller FS correction compared to that of a cubic system. 
This is discussed further in the next Section. 

These applications demonstrate that the new FS error 
correction scheme can consistently remove most of the FS 
error and accelerate the size convergence of MB results 
to infinite size limit. 



V. DISCUSSION 

In this section we discuss several aspects of the new FS 
correction method, the parametrized LSDA functional, 
and possible further improvements. 

In Sec. [Hi] the FS LSDA XC functional was separated 
into exchange £ x (r s , L, £) [Eqs. (|ll[) - (fT5)) ] and correlation 
e c (r s ,L, C) [Eqs. ([TBI) arid ([fT])] contributions. This facil- 
itates the modification of existing DFT computer codes 
to generate the FS corrections. While useful, this sepa- 
ration somewhat complicates the FS treatments, which 
are based on XC functionals derived from the HEG. 

For example, while the total FS HEG XC energy scales 
as 1/L 3 , e x (r s ,L, £) scales as 1/L 2 , since the Hartree- 
Fock spectrum is gapless in this case. The separation 
thus induces a canceling 1/L 2 term in e c (r s , L, (). Partly 
for this reason, we adopted the PZ spin interpolation for- 
mula given by Eq. (fTTj) for e x (r s , L,(), rather than the 
exact relation in Eq. (|10[) , since it makes the treatment 



of exchange consistent with that of the correlation, for 
which there is no analogy of Eq. (|10j) . and ensures can- 
cellation of the 1/L 2 dependence. [As mentioned, in the 
limit L — » oo, Eq. (fTTj) reduces to the exact relation in 
Eq. (flO|) .] Since the FS exchange energy for spin un- 
polarized system is discontinuous,— another reason for 
using Eq. (jllj) is that it avoids multiple discontinuities 
in a partially polarized system, as discussed in Sec. 1111 Al 
The continuity density j x in Eq. (fT3j) was chosen for sim- 
plicity. It could be made polarization-dependent, but the 
expressions would be more complex and not as easy to 
implement into DFT codes. As discussed in Sec. IIII CI 
FS corrections are not sensitive to the precise value of 
jx, especially for extended systems. 

Similar considerations pertain to the treatment of cor- 
relation. Since e c (r s , L,£) is continuous, it was easy to 
use different values for 7 C (0) and 7 C (1) in Eq. (|T6|) . as 
described in Sec. IIII CI when the PZ interpolation form 
is used. If the PW interpolation of Eq. (fl8|) were to 
be used, cancellation of the 1/L 2 term would require 
more complicated expressions. Interpolation formulas, 
such as that of PZ, allow us to map spin unpolarizcd 
and fully-polarized results onto the correlation energy at 
an arbitrary polarization. For the fully polarized HEG, 
the available finite-size data from QMC were not as de- 
tailed as for the unpolarized case. Our fits of the fully- 
polarized e c (r s , L,p = 1) correlation functional arc based 
on AFQMC calculations of finite systems, mostly for r s 
greater than about 1.5 - 2. While our correlation en- 
ergy fits are acceptable (see Figs. [2] and [3]), more QMC 
calculations for small r s might lead to improvements. 

Extensions of our FS error correction scheme are 
possible, for example based on orbital-dependent den- 
sity functionals,— such as hybrid DFT functional cal- 
culations. Hybrid functional s 42 ' 44 augment the DFT 
exchange-correlation energy with an exact exchange term 
calculated from HF theory. Hybrid functionals have be- 
come one of the most promising approaches for overcom- 
ing some of the shortcomings of local or scmilocal DFT 
approximations." 4 ^ Some complications would need to be 
addressed. Popular hybrid functionals such as PBEO 4 ^ 
use a fraction ~ 20 % of exact exchange and DFT ex- 
change for the remainder. In applications to semicon- 
ducting and insulating materials, the HF exchange con- 
verges as 1/L 3 with PBC, since the single-particle spec- 
trum has a gap. [The leading behavior is 1 /L 3 after a 1/L 
Madclung-like dependence is removed ) 45 ' 46 as is done in 
QMC calculations.—] Since DFT FS exchange scales as 
1 /L 2 for the HEG, the fraction of the 1 / L 2 term in the FS 
correlation functional would have to be adjusted in order 
to obtain a net 1/L 3 behavior. The resulting e c (r s , L, £) 
would lead to a non-continuous behavior with increasing 
r s , which might require further modifications of the fits. 

Finally, our FS XC functional is based on HEG calcu- 
lations using cubic supercells of edge L. For non-cubic 
applications, L is determined by the volume of the su- 
percell. In the HEG, spatially isotropic supercells are 
favored. For example, a small cubic to tetragonal distor- 
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tion will cause the total energy to become more positive, 
scaling as ~ (1 — C * r/ 2 ), where C > is a constant and 
r\ describes the aspect ratio c/a=l + r;of the tetrago- 
nal supercell. Therefore our FS XC functional tends to 
provide an upper limit to the magnitude of the FS correc- 
tion for non-cubic supercells. This may have contributed 
to the over-correction seen in Fig. [71 where anisotropic 
MnO and FeO supercells were used. 

VI. SUMMARY 

Many-body calculational methods such as quantum 
Monte Carlo can potentially provide accurate results for 
extended systems where electron interaction effects arc 
important. These calculations are computationally more 
costly than independent-electron calculations, and are 
always accompanied by more significant FS errors. In 
this work, we have provided a framework within DFT to 
understand and estimate the FS errors in MB calcula- 
tions of spin-polarized systems. A FS LSDA functional 
is parametrized which leads to a simple, post-processing 
correction method. The method is designed to approxi- 
mately include 2B FS corrections in FS DFT calculations 



of systems with arbitrary spin polarization. The correc- 
tions to total energy, dissociation energy, and ionization 
energy calculations of several typical systems show a sig- 
nificant improvement on the convergence. The formalism 
of the method and the parametrization of the FS func- 
tional are presented in detail. The strengths and weak- 
nesses of the approach are discussed. There are several 
directions in which the method can be further developed. 
The method can be easily incorporated in any standard 
DFT computer program for periodic systems, and can be 
used to correct any MB results. It is expected that the 
method will find applications in the study of a variety of 
realistic systems with magnetic ordering. 
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